Clear["Global`*"]
(* This command wipes everything out from the current session of Mathematica. *)

CFT and gravity, solutions week 6

OPE scalar case (questions 1c. and 1d.)

Scalar product and differentiation (copied from week 3)

Define (Subscript[g, μν]), a scalar product · and rules for contractions

Attributes[CenterDot]={Orderless};
(a_+b_)·c_:=a·c+b·c
a_·(α_ b_):=α (a·b)/;NumberQ[α]
Attributes[g]=Orderless;
dot:=((Expand[#]//. {
z_[μ_] y_[μ_]:>z·y,((z_[μ_])^2):>z·z,g[μ_,ν_]z_[ν_]:>z[μ],g[μ_,α_]g[α_,ν_]:> g[μ,ν] ,((g[μ_,ν_])^2):>d,g[μ_,μ_]:>d})//Simplify)&

Warning: it does not handle mistakes in the repeated indices: see first example below. So please contract indices properly in the following.

x[μ]x[α] y[μ]x[μ]//dot
x[μ]x[α] y[μ]x[α]//dot
x[μ](x[μ]+y[μ])//dot
x[μ](x[ν]+y[ν])g[μ,ν]//dot
g[μ,ν]g[ν,μ]//dot

Define (∂_μ)

der[x_][μ_][expr_] takes three arguments: the point x, the direction μ and the expression to be derived. We use the command SeriesCoefficient that knows already about linearity and Leibnitz rule:

der[x_][μ_][A_]:=(A/.{
x[α_]:> x[α]+ϵ g[μ,α],x·x:> x·x+2ϵ x[μ],x·y_:> x·y+ϵ y[μ]})//SeriesCoefficient[#,{ϵ,0,1}]&

der[x][μ][x·x]
der[x][μ][x·y]
der[x][μ][f[x·y]]

Constructing the basic ingredients involved in the equation that fixes the OPE

ClearAll[Fder]

(*simplify[X_]:=Simplify[ (Expand[X]/.{w·w->1})]*)
d1[X_]:=(der[y][μ][X]x[μ]//dot)/.{w·w->1}
d2[X_]:=(der[y][ν][der[y][μ][X]]g[μ,ν]x·x//dot)/.{w·w->1}
Fder[0,0]=(1-2y·w+y·y)^(-Δ3);
Fder[n_,m_]:=(Fder[n,m]=d1[Fder[n-1,m]])/;n>0;
Fder[n_,m_]:=(Fder[n,m]=d2[Fder[n,m-1]])/;m>0;

Notice the structure f[x_]:=f[x]=stuff. This is called caching. The effect of the second equal sign is that when a function is called with argument x0, the value it returns is stored in the variable f[x0]. This way if you call again f[x0], mathematica does not compute it anymore, it just returns the value previously stored. You can read more about this looking for the following link on the help: tutorial/FunctionsThatRememberValuesTheyHaveFound

Check

Fder[2,3]-Nest[d1,Nest[d2,(1-2y·w+y·y)^(-Δ3),3],2]

Determining the OPE coefficients of descendants (question 1c.)

M=10;
RHS=(1-2x·w+x·x)^((Δ2-Δ1-Δ3)/2);
LHS=Sum[a[n,m]Fder[n,m],{m,0,M/2},{n,0,M-2m}]/.{a_·y->0}//Simplify;
vars=Table[a[n,m] ,{m,0,M/2},{n,0,M-2m}]//Flatten;
LHS-RHS/.{x·x-> t^2,x·w-> t b};
Series[%,{t,0,M}]//Normal;
CoefficientList[%,{t,b}]//Flatten;
Timing[data=Solve[%==0,vars]//Factor;]

Guessing the general term (question 1d.)

FindSequenceFunction is a good way to start.

Table[a[n,0]/.data[[1]],{n,0,M}];
FindSequenceFunction[%,n+1]//Simplify

But it may take too much time, so let’s try to guess.

Table[a[0,m] /(Pochhammer[(1/2) (Δ1-Δ2+Δ3),m]Pochhammer[(1/2)(-Δ1+Δ2+Δ3),m])/.data[[1]],{m,0,M/2}]//Simplify
Table[a[0,m]  Pochhammer[Δ3,2m] /(Pochhammer[(1/2) (Δ1-Δ2+Δ3),m]Pochhammer[(1/2)(-Δ1+Δ2+Δ3),m])/.data[[1]],{m,0,M/2}]//Simplify
Table[a[0,m]  Pochhammer[(2Δ3-d+2)/2,m]Pochhammer[Δ3,2m] /(Pochhammer[(1/2) (Δ1-Δ2+Δ3),m]Pochhammer[(1/2)(-Δ1+Δ2+Δ3),m])/.data[[1]],{m,0,M/2}]//Simplify
Table[a[0,m]Pochhammer[(2Δ3-d+2)/2,m]Pochhammer[Δ3,2m] /(Pochhammer[(1/2) (Δ1-Δ2+Δ3),m]Pochhammer[(1/2)(-Δ1+Δ2+Δ3),m])/.data[[1]],{m,0,M/2}]//Simplify
FindSequenceFunction[1/%,m+1]//FunctionExpand
Table[a[0,m](((-4))^m) m!Pochhammer[(2Δ3-d+2)/2,m]Pochhammer[Δ3,2m] /(Pochhammer[(1/2) (Δ1-Δ2+Δ3),m]Pochhammer[(1/2)(-Δ1+Δ2+Δ3),m])/.data[[1]],{m,0,M/2}]//Simplify

Then we put together, with a bit of trial and error, and get

aguess[n_,m_]:=Pochhammer[((-Δ1+Δ2+Δ3)/(2 )),m]/((-1)^(m) n!4^m m! Pochhammer[Δ3,n+2m]Pochhammer[(2Δ3-d+2)/2,m])Pochhammer[((Δ1-Δ2+Δ3)/(2 )),n+m]
Table[a[n,m]/aguess[n,m] ,{m,0,M/2},{n,0,M-2m}]/.data[[1]]//Factor